knitr::opts_chunk$set(echo = TRUE)
list.of.packages <- c("gsheet", "tidyverse", "janitor", "plotly", "glmmTMB", "metafor", "broom.mixed", "car") #add new libraries here
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)
# Load all libraries
lapply(list.of.packages, FUN = function(X) {
do.call("require", list(X))
})
sessionInfo()
NOTE: data is read directly from the GoogleSheet using a share link that was set to “anyone with a link can view”
# Read in data from GoogleSheet
data.fert <- as_tibble(gsheet2tbl('https://docs.google.com/spreadsheets/d/111SuH548Et6HDckjbQjtYk-B8A5VfNkUYZgcUNRPxxY/edit?usp=sharing'))
## Warning: Missing column names filled in: 'X22' [22]
#replace NR (not reported) with NA and convert columns to factor / numeric where needed
data.fert <- data.fert %>%
na_if("NR") %>%
mutate_at(c('Phylum', 'Study', 'Taxonomic Group', 'Common name', 'Latin name', 'Error statistic'), as.factor) %>%
mutate_at(c('pH Experimental', 'pH Control', 'pCO2 Experimental', 'pCO2 Control', 'Ave. Fert. % @ pH', 'Error % @ pH', '# Trials @ pH', 'Insemination time'), as.numeric) %>%
clean_names() # fill in spaces with underscores for column names
## Warning: NAs introduced by coercion
## Warning: NAs introduced by coercion
#data.fert %>% select(study, ave_fert_percent_p_h, error_percent_p_h, error_statistic, number_trials_p_h) %>% View()
ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) +
geom_point(size=1.5, width=0.02) +
facet_wrap(~phylum, scale="free") +
geom_smooth(method="lm", se=TRUE, aes(fill=taxonomic_group)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum, linear") +
theme_minimal()
)
## Warning: Ignoring unknown parameters: width
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
#ggplotly(
data.fert %>%
ggplot(mapping=aes(x=p_h_experimental, y=ave_fert_percent_p_h, group=taxonomic_group, col=taxonomic_group, text=`common_name`)) +
geom_point(size=1, width=0.02) +
facet_wrap(~phylum, scale="free") +
geom_smooth(method="lm", se=TRUE, formula=y ~ poly(x, 2, raw=TRUE), aes(fill=taxonomic_group)) +
ggtitle("Fertilization Rate ~ pH exposure by phylum, polynomial") +
theme_minimal()
## Warning: Ignoring unknown parameters: width
## Warning: Removed 25 rows containing non-finite values (stat_smooth).
## Warning in predict.lm(model, newdata = new_data_frame(list(x = xseq)), se.fit =
## se, : prediction from a rank-deficient fit may be misleading
## Warning: Removed 25 rows containing missing values (geom_point).
#)
Upper 95% CI = Mean + 1.96*SE; I recorded the difference between the upper 95%CI and the mean (Upper 95%CI - Mean). To convert I will use:
SE = (Upper 95%CI - Mean) / 1.96
SD = ((Upper 95%CI - Mean) / 1.96) * sqrt(n)
SE= SD/sqrt(n); where n=sample size. To convert I will use that equation (since I recorded sample size, i.e. number of trials at each pH)
SD = SE*sqrt(n)
I will use that statistic as-is (no conversion), thereby assuming it is SE.
data.fert <- data.fert %>%
mutate(SE = case_when(error_statistic == "SD" ~ error_percent_p_h/sqrt(number_trials_p_h),
error_statistic == "95% CI" ~ error_percent_p_h/1.96,
is.na(error_statistic) ~ error_percent_p_h,
error_statistic == "SE" ~ error_percent_p_h))
data.fert <- data.fert %>%
mutate(SD = case_when(error_statistic == "SE" ~ error_percent_p_h*sqrt(number_trials_p_h),
error_statistic == "95% CI" ~ (error_percent_p_h/1.96)*sqrt(number_trials_p_h),
is.na(error_statistic) ~ error_percent_p_h,
error_statistic == "SD" ~ error_percent_p_h))
data.fert <- data.fert %>%
mutate(pH_delta = p_h_control-p_h_experimental)
data.fert <- data.fert %>%
mutate(ave_fert_proport = case_when(ave_fert_percent_p_h <= 100 ~ ave_fert_percent_p_h/100,
ave_fert_percent_p_h > 100 ~ 1))
Transormation equation source: https://cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf, " … if y also assumes the extremes 0 and 1, a useful transformation in practice is (y · (n − 1) + 0.5)/n where n is the sample size (Smithson and Verkuilen 2006)."
issue is that a few studies only had 1 trial per pH, so the transformation results in NA values
# data.fert$ave_fert_proport.t <- data.fert$ave_fert_proport*((data.fert$number_trials_p_h-1) + 0.5) / data.fert$number_trials_p_h
data.fert$ave_fert_proport.t <- data.fert$ave_fert_proport - 0.001
data.fert$ave_fert_proport.t %>% hist()
ggplot(data.fert, aes(group=phylum, col=phylum)) + geom_density(aes(ave_fert_proport.t))
## Warning: Removed 13 rows containing non-finite values (stat_density).
# How many studies per ~phylum?
data.fert %>%
select(phylum, study) %>%
distinct(phylum, study) %>%
group_by(phylum) %>% count()
## # A tibble: 4 x 2
## # Groups: phylum [4]
## phylum n
## <fct> <int>
## 1 Cnidaria 4
## 2 Crustacea 5
## 3 Echinodermata 12
## 4 Mollusca 18
# How many studies per taxonomic group?
data.fert %>%
select(phylum, study, taxonomic_group) %>%
distinct(phylum, study, taxonomic_group) %>%
group_by(taxonomic_group) %>% count()
## # A tibble: 10 x 2
## # Groups: taxonomic_group [10]
## taxonomic_group n
## <fct> <int>
## 1 abalone 2
## 2 clam 5
## 3 copepod 3
## 4 coral 4
## 5 mussel 4
## 6 non-copepod 2
## 7 non-urchin 2
## 8 oyster 5
## 9 scallop 3
## 10 sea urchin 11
weights <- metafor::escalc(measure='MN',
mi=data.fert$ave_fert_percent_p_h,
sdi = data.fert$SD,
ni=data.fert$number_trials_p_h, options(na.action="na.pass"))
data.fert$w <-weights$vi
test1 <- glmmTMB(ave_fert_proport ~ p_h_experimental + taxonomic_group/phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
test2 <- glmmTMB(ave_fert_proport ~ p_h_experimental*phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test3 <- glmmTMB(ave_fert_proport ~ p_h_experimental:phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test4 <- glmmTMB(ave_fert_proport ~ p_h_experimental + phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test5 <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test6 <- glmmTMB(ave_fert_proport ~ 1 + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test7 <- glmmTMB(ave_fert_proport ~ phylum + (1|study), data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
test8 <- glmmTMB(ave_fert_proport ~ p_h_experimental, data=drop_na(data.fert, w), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
AIC(test1, test2, test3, test4, test5, test6, test7, test8) #test2 smallest AIC.
## Warning in AIC.default(test1, test2, test3, test4, test5, test6, test7, : models
## are not all fitted to the same number of observations
## df AIC
## test1 42 NA
## test2 9 1169.716
## test3 6 1208.641
## test4 6 1206.182
## test5 3 1210.887
## test6 2 1450.084
## test7 5 1440.438
## test8 2 3511.136
#test differene between models
anova(test3, test2)
## Data: drop_na(data.fert, w)
## Models:
## test3: ave_fert_proport ~ p_h_experimental:phylum + (1 | study), zi=~0, disp=~1
## test2: ave_fert_proport ~ p_h_experimental * phylum + (1 | study), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## test3 6 1208.6 1230.7 -598.32 1196.6
## test2 9 1169.7 1202.9 -575.86 1151.7 44.924 3 9.601e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Examine model test2
car::Anova(test2) #phylum, pH, & phylum:pH sign. factors
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 156.796 1 < 2.2e-16 ***
## phylum 10.563 3 0.01434 *
## p_h_experimental:phylum 42.497 3 3.147e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(test2)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental * phylum + (1 | study)
## Data: drop_na(data.fert, w)
## Weights: w
##
## AIC BIC logLik deviance df.resid
## 1169.7 1202.9 -575.9 1151.7 285
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.9459 0.9726
## Number of obs: 294, groups: study, 31
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.1350 5.3425 1.148 0.2508
## p_h_experimental -0.7656 0.6731 -1.137 0.2554
## phylumCrustacea -12.5491 5.5114 -2.277 0.0228 *
## phylumEchinodermata -22.6771 5.5012 -4.122 3.75e-05 ***
## phylumMollusca -13.7549 6.0574 -2.271 0.0232 *
## p_h_experimental:phylumCrustacea 1.6380 0.6958 2.354 0.0186 *
## p_h_experimental:phylumEchinodermata 2.9631 0.6942 4.268 1.97e-05 ***
## p_h_experimental:phylumMollusca 1.8986 0.7692 2.468 0.0136 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Generate estimates & confidence intervals (log likelihood)
confint(test2)
## 2.5 % 97.5 % Estimate
## cond.(Intercept) -4.3361108 16.6060479 6.1349685
## cond.p_h_experimental -2.0847698 0.5536213 -0.7655743
## cond.phylumCrustacea -23.3512502 -1.7470097 -12.5491300
## cond.phylumEchinodermata -33.4593535 -11.8949260 -22.6771397
## cond.phylumMollusca -25.6272057 -1.8825189 -13.7548623
## cond.p_h_experimental:phylumCrustacea 0.2741953 3.0017240 1.6379596
## cond.p_h_experimental:phylumEchinodermata 1.6025270 4.3237081 2.9631176
## cond.p_h_experimental:phylumMollusca 0.3910182 3.4062259 1.8986221
## study.cond.Std.Dev.(Intercept) 0.7333418 1.2899079 0.9725962
# Instpect residuals ~ fitted values
aa5 <- augment(test2, data=drop_na(data.fert, w))
## Warning in mr - fitted(object): longer object length is not a multiple of
## shorter object length
#ggplotly(
ggplot(aa5, aes(x=.fitted,y=.resid)) +
geom_point() +
geom_smooth()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 12 rows containing non-finite values (stat_smooth).
## Warning: Removed 12 rows containing missing values (geom_point).
#)
ph.min.max <- drop_na(data.fert, w) %>%
select(phylum, p_h_experimental) %>%
group_by(phylum) %>%
summarize(min=min(p_h_experimental, na.rm=TRUE), max=max(p_h_experimental, na.rm=TRUE))
phylum.list <- list()
for (i in 1:nrow(ph.min.max)) {
phylum.list[[i]] <- data.frame(ph=c(seq(from=as.numeric(ph.min.max[i,"min"]),
to=as.numeric(ph.min.max[i,"max"]),
by=0.01)),
phylum=rep(c(ph.min.max[i,"phylum"])))
}
new.data <- bind_rows(phylum.list) %>% purrr::set_names(c("p_h_experimental", "phylum"))
new.data$study <- NA
new.data$w <- NA
predict.test.df <- predict(test2, newdata = new.data, se.fit = TRUE, type="response")
predict.test.df.df <- predict.test.df %>%
as.data.frame() %>%
cbind(new.data)
#scales::show_col(c("#e41a1c","#4daf4a","#ff7f00","#984ea3",'#377eb8'))
# Data with beta regression model fit
#ggplotly(
ggplot() +
geom_jitter(data=drop_na(data.fert, w), aes(x=p_h_experimental, y=ave_fert_proport, group=phylum, col=phylum), size=1.2, width=0.03) +
facet_wrap(~phylum, scales="free") + theme_minimal() +
ggtitle("% fertilization ~ pH with binomial-regression model predictions") +
xlab("Experimental pH") + ylab("Fertilization %") +
scale_color_manual(name=NULL, values=c("#e41a1c","#ff7f00","#4daf4a",'#377eb8')) +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.test.df.df, aes(x=p_h_experimental, y=fit), col="gray50") + #, col=phylum
geom_ribbon(data = predict.test.df.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit, fill=phylum), linetype=2, alpha=0.1, col="gray50") #)
## Warning: Removed 12 rows containing missing values (geom_point).
Fertilization success (%) by experimental pH across marine taxa examined in this review. Meta-analysis was performed using a binomial regression model, and indicates that fertilization success decreases with pH across Cnidaria (4 studies), Crustacea (5 studies), Echinodermata (12 studies), and Mollusca (18 studies). Each point reflects the average % fertilization reported by one study at an experimental pH.
model.cnidaria <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.cnidaria)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 0.7736 1 0.3791
summary(model.cnidaria)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Cnidaria")
## Weights: w
##
## AIC BIC logLik deviance df.resid
## 155.8 161.5 -74.9 149.8 46
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.009966 0.09983
## Number of obs: 49, groups: study, 3
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.4307 5.2035 0.852 0.394
## p_h_experimental -0.5806 0.6601 -0.880 0.379
predict.cnidaria <- predict(model.cnidaria, newdata = subset(new.data, phylum=="Cnidaria"), se.fit = TRUE, type="response")
predict.cnidaria.df <- predict.cnidaria %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Cnidaria"))
ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>% filter(phylum=="Cnidaria"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, color="#e41a1c") +
ggtitle("Cnidaria") +
xlab("Experimental pH") + ylab("Fertilization %") +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.cnidaria.df, aes(x=p_h_experimental, y=fit), col="#e41a1c") +
geom_ribbon(data = predict.cnidaria.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#e41a1c")
model.crustacea <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), binomial(link = "logit"), na.action=na.exclude, weights=w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.crustacea)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 26.088 1 3.263e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.crustacea)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Crustacea")
## Weights: w
##
## AIC BIC logLik deviance df.resid
## 119.6 123.2 -56.8 113.6 21
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 1.837 1.355
## Number of obs: 24, groups: study, 5
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.0292 1.5082 -4.661 3.15e-06 ***
## p_h_experimental 0.9159 0.1793 5.108 3.26e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.crustacea <- predict(model.crustacea, newdata = subset(new.data, phylum=="Crustacea"), se.fit = TRUE, type="response")
predict.crustacea.df <- predict.crustacea %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Crustacea"))
ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>% filter(phylum=="Crustacea"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#ff7f00") +
ggtitle("Crustacea") +
xlab("Experimental pH") + ylab("Fertilization %") +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.crustacea.df, aes(x=p_h_experimental, y=fit), col="#ff7f00") +
geom_ribbon(data = predict.crustacea.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#ff7f00")
model.echinodermata <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), binomial(link = "logit"), na.action=na.exclude, weights = w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.echinodermata)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 154.04 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.echinodermata)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Echinodermata")
## Weights: w
##
## AIC BIC logLik deviance df.resid
## 700.2 709.2 -347.1 694.2 142
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.9716 0.9857
## Number of obs: 145, groups: study, 11
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -15.8468 1.3623 -11.63 <2e-16 ***
## p_h_experimental 2.1430 0.1727 12.41 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.echinodermata <- predict(model.echinodermata, newdata = subset(new.data, phylum=="Echinodermata"), se.fit = TRUE, type="response")
predict.echinodermata.df <- predict.echinodermata %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Echinodermata"))
ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>% filter(phylum=="Echinodermata"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#4daf4a") +
ggtitle("Echinodermata") +
xlab("Experimental pH") + ylab("Fertilization %") +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.echinodermata.df, aes(x=p_h_experimental, y=fit), col="#4daf4a") +
geom_ribbon(data = predict.echinodermata.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#4daf4a")
model.mollusca <- glmmTMB(ave_fert_proport ~ p_h_experimental + (1|study), data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), binomial(link = "logit"), na.action=na.exclude, weights=w)
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
Anova(model.mollusca)
## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: ave_fert_proport
## Chisq Df Pr(>Chisq)
## p_h_experimental 9.4414 1 0.002121 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(model.mollusca)
## Family: binomial ( logit )
## Formula: ave_fert_proport ~ p_h_experimental + (1 | study)
## Data: data.fert %>% drop_na(w) %>% filter(phylum == "Mollusca")
## Weights: w
##
## AIC BIC logLik deviance df.resid
## 198.1 205.1 -96.1 192.1 73
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## study (Intercept) 0.818 0.9045
## Number of obs: 76, groups: study, 14
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.563 2.811 -2.690 0.00714 **
## p_h_experimental 1.128 0.367 3.073 0.00212 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
predict.mollusca <- predict(model.mollusca, newdata = subset(new.data, phylum=="Mollusca"), se.fit = TRUE, type="response")
predict.mollusca.df <- predict.mollusca %>%
as.data.frame() %>%
cbind(subset(new.data, phylum=="Mollusca"))
ggplot() +
geom_jitter(data=data.fert %>% drop_na(w) %>% filter(phylum=="Mollusca"), aes(x=p_h_experimental, y=ave_fert_proport), size=1.2, width=0.03, col="#377eb8") +
ggtitle("Mollusca") +
xlab("Experimental pH") + ylab("Fertilization %") +
theme_minimal() +
coord_cartesian(ylim = c(0, 1), xlim =c(6,8.5)) +
geom_line(data = predict.mollusca.df, aes(x=p_h_experimental, y=fit), col="#377eb8") +
geom_ribbon(data = predict.mollusca.df, aes(x=p_h_experimental, ymin=fit-se.fit, ymax=fit+se.fit), linetype=2, alpha=0.1, col="#377eb8")
## Warning: Removed 12 rows containing missing values (geom_point).